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Dynamic Smagorinsky model on anisotropic grids 

By A. Scotti 1 , C. Meneveau 1 AND M. Fatica 2 


Large Eddy Simulation (LES) of complex-geometry flows often involves highly 
anisotropic meshes. To examine the performance of the dynamic Smagorinsky 
model in a controlled fashion on such grids, simulations of forced isotropic turbu- 
lence are performed using highly anisotropic discretizations. The resulting model 
coefficients are compared with a theoretical prediction (Scotti ct al , 1993). Two 
extreme cases are considered: pancake-like grids, for which two directions are poorly 
resolved compared to the third, and pencil-like grids, where one direction is poorly 
resolved when compared to the other two. For pancake-like grids the dynamic model 
yields the results expected from the theory (increasing coefficient with increasing as- 
pect ratio), whereas for pencil-like grids the dynamic model does not agree with the 
theoretical prediction (with detrimental effects only on smallest resolved scales). A 
possible explanation of the departure is attempted, and it is shown that the problem 
may be circumvented by using an isotropic test-filter at larger scales. 

Overall, all models considered give good large-scale results, confirming the gen- 
eral robustness of the dynamic and eddy-viscosity models. But in all cases, the 
predictions were poor for scales smaller than that of the worst resolved direction. 


1. Introduction 

Since its introduction in the 196CTs, a goal of LES has been to simulate complex 
turbulent flows. A complex flow is, by definition, characterized by regions were the 
physics of turbulence change, e.g., from homogeneous turbulence far from bound- 
aries to near wall turbulence, eic. To capture the full gamut with a simple subgrid 
model without having to adjust constants in an ad hoc manner every time was a 
serious problem until recently. The introduction of the dynamic model (Germano et 
a/., 1991) to dynamically calculate the parameter(s) of the modeled sub-grid stress 
was a significant step towards making LES of complex flows possible without ad hoc 
adjustments. This model is able to self-adjust to the large scale flow in the correct 
fashion, for instance, shutting itself down near walls or in regions where the flow 
relaminarizes. 

As a result, it has become possible to apply LES to study flows of increasing 
complexity (e.g. Akselvoll and Moin 1996 or see in this same volume Chan and 
Mittal, and Haworth and Jansen), which in turn requires the use of complex grids, 
either structured or unstructured. Complicated grid geometries in conjunction with 
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the dynamic model raise several questions. Consider, as an example, the flow past 
a 3-D bluff body: near the object, one needs to refine the grid in the spanwise direc- 
tions. For a structured mesh, far downstream, the grid may be greatly expanded in 
the streainwise direction. Therefore, the grid can be strongly anisotropic, with the 
elements of the grid looking like sheets or pencils, depending on the kind of refine- 
ment imposed upstream. Hence, in the far- wake region one may have a situation 
where the turbulence is nearly isotropic, whereas the computational grid is highly 
anisotropic. 

In LES, the grid filter is dictated by the computational mesh used to solve the 
equations (although, for methods other than spectral, it is difficult to give a precise 
definition of the filtering operator associated with a given discretization). Since 
classical eddy-viscosity models need as input a length-scale which is usually associ- 
ated with the scale at which the filter operates, the problem arises in defining this 
length when, as a result of the anisotropy of the grid, the filter is defined by more 
than one length scale. For the Srnagorinsky model, this problem was considered 
first by Deardoff (1970) and later by Schumann (1975), Lilly (1988) and Scotti et 
al. (1993), although the last two papers were only theoretical treatments. 

On the other hand, other models such as the dynamic model do not in principle 
require a length scale to be specified. The question then arises whether the dynamic 
model is able to correctly simulate isotropic turbulence on anisotropic grids. The 
main goal of this work is to examine this question. 

This issue is also of theoretical interest since, from the point of view of interaction 
among modes, local triadic interactions at small scales are fully available only to a 
limited amount of modes. Thus the small scales are exposed to a dynamic which 
is not the one typical of 3-D turbulence. It is natural then to expect that the SGS 
stress tensor should incorporate a correction originating only from the anisotropy 
of the grid. 

The paper is organized as follows: in section 2 we briefly summarize the main 
result of Scotti et al. (1993) and set the notation that will be used throughout 
the paper; in Section 3 we discuss the simulations and how the results of different 
models will be compared. In showing the results, we have considered two categories 
of grids: pancake- like, when one direction is much better resolved than the other 
two, and pencil- like, when two directions are much better resolved than the third. 
Section 4 presents the results. Finally, in Section 5 a summary and discussion of 
the results is given. 

2. Srnagorinsky model on anisotropic grids 

In this section, the results of Scotti et al. (1993) are briefly recalled. They are 
based on the assumption that the turbulence is isotropic and homogeneous, and 
that the largest and smallest scales at which the filter operates still lie within the 
inertial range. One begins by writing the Smagorinsky model as 

m = - 2[£( A j , A 2 , A 3 )] 2 [25f m ] 1 ,2 Si] . (1) 

Here Ai, A 2 and A 3 are the dimensions of the computational cell. For notational 
convenience and without lack of generality, let us assume Ai < A 2 < A 3 . The 
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equivalent filter, via a collocation rule, is assumed to be a sharp cut-off filter in 
Fourier space, which corresponds to setting to zero all the modes outside the region 
B = {\ki \ < tt/A 1? 1*2 1 < 71 -/A 2 , | A: 3 1 < 7 r / A 3 } , leaving the others unmodified. 

By invoking an argument used first by Lilly (1967) an expression for L(Ai, A 2 , A 3 ) 
was derived by requiring that 


£ — < T{j S%j >, 

replacing r tJ with the model and computing moments of the strain-rate tensor, 
assuming that the velocity field is characterized by a Kolmogorov isotropic spectrum 
on all resolved modes. 

Introducing A eq = (A 1 A 2 A 3 ) 1 / 3 , L(Ai, A 2 , A 3 ) can then be written as 


L(Ai,A 2,A*) = C , <9 A e9 /(ai,a 2 ), (2) 

where a\ = Ai / A 3 and a 2 = A 2 /A 3 are the two aspect ratios of the grid, and f > 1 
is a function equal to one if both ratios are equal to unity. C s is the traditional 
Smagorinsky coefficient, which depends on the value of the Kolmogorov constant. 

After evaluating the function /, a compact approximation for the result was given 
by Scotti et ai (1993) 

/(a 1 , « 2 ) — cosh v / 4/27((loga 1 ) 2 - log a! log 02 + (loga 2 ) 2 ). (3) 


Incidentally, we remark that the fact that / ~ 1 for aspect ratios close to unity 
justifies the practice introduced by Deardoff (1970) of using A eg as length scale, 
at least for aspect ratios close to unity. In the dynamic version of this model, 
with grid filtering denoted by tilde and test filtering by an overbar, the length- scale 
X(Ai,A 2 , A 3 ) is computed according to 


where 


2[I(A,,A 2 ,A 3 )] 2 


< LjjMjj > 

< MijMij > ’ 


Lij = Ujiii — UiUj, 


( 4 ) 


(4a) 


and 


M tJ = 



A e i 1 /(Ql,Q2) 

A e , /(di,a 2 ) 


2 


nsLrs 



(46) 


where we have made use of Eq. (2). If both test and grid filter have the same aspect 
ratios then Eq. (4) is closed; otherwise we can use Eq. (3) to compute / and check 
a posteriori its consistency. 
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3. Approach and validation 

We run LES of isotropic turbulence in a box of side 2ir with periodic boundary 
conditions. Turbulence is maintained by a forcing f that forces the largest modes 
(Jt < 2) with an intensity such that the energy injection rate f • u is fixed at a 
constant value £ = 1.0. The numerical scheme is the same as in Vincent and 
Meneguzzi (1991) and Briscolini and Santangelo (1994). It uses Adam-Bashforth 
2 for time advancing, with At — 0.001. The nonlinear terms, written in rotational 
form, are evaluated pseudospectrally. Appendix A examines dealiasing for the AB2 
scheme. The grids have mesh sides (Ai,A 2 ,A 3 ), with A 3 > max{Ai, A 2 }, and 
aspect ratios a\ = Aj/A 3i a 2 = A 2 /A 3 ranging from 1 to 1/16. Grid filtering 
was performed with a sharp spectral cut-off setting to zero the modes outside the 
ellipsoid B = {k G R 3 | (AqAi ) 2 + (k 2 A 2 ) 2 + (A 3 A 3) 2 < 8/9tt 2 }, which has the 
advantage of partially removing aliasing errors (see appendix A). Test filtering was 
done at a scale twice as large in all directions. 

For comparison, computations were performed using the classical non-dynamic 
Smagorinsky model with the Deardoff length scale and C 2 = 0.026, as well as with 
the Smagorinsky model corrected after Scotti at al. (1993) including f(aj,a 2 ) as 
evaluated from Eq. (3). In all cases the initial condition is assumed to be a random 
Gaussian field with A ~ 4 5 / 3 spectrum, random phase, and total kinetic energy equal 
to unity. 

We wish to compare both large scale properties, such as total kinetic energy, 
derivative skewness in the worst resolved direction, and small scale properties, such 
as energy spectra near cut-off scale and the skewness in the best resolved direction 
(which is sensitive to the details of the small scales). 

For isotropic turbulence we know that the spectral tensor in the inertial range is 
given bv 

Q li (k)=<ti < (k)« > (-k)>=(47r)- 1 CK£ 2 / 3 fc- 11 / 3 P 1 J (k), (5) 

where e is the average dissipation, Ck is the Kolmogorov constant, and P,j(k) is 
the projector on the space orthogonal to k. Also, we know that the skewness of the 
derivative is 0(— . 5 ), although for LES the value attained is typically smaller due 
to the incomplete resolution of the small scales. We will compute the skewness in 
the a— direction, defined as S a =< ( du Q /dx a ) 3 > / < ( du Q /dx a ) 2 > 3 / 2 . 

Due to the anisotropy of the grid, it is better to study 1 -D premultiplied spectra, 
defined as 

, J B 2n S - 2 /*k n / 3 Q„(k)dk 2 dk 3 
(l) ~ f B dk 2 dk 3 

For ideal Kolmogorov turbulence, where the spectral tensor is given by Eq. (5), 
C(A*i) is a constant equal to the Kolmogorov constant Ck ~ 1.6. 

4. Results 

To obtain a self-consistent estimate for the Smagorinsky constant C 3 , we first run 
LES with the dynamic model with isotropic spherical test and grid filter on a 32 3 
grid. After an initial transient the value stabilizes at C 2 = 0.023 ± 5%. Next, we 
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FIGURE 1. (a) Time traces of /dyn(ai,a 2 ) as generated by the dynamic model 

during LES of forced turbulence on anisotropic grids. : aspect ratios a\ = 

a 2 = 1/8; : ai = a 2 = 1; : ai = 1/16, a 2 = 1. (b) values of time 

averages of /dyn(«i,«2) computed between 400 < t < 800, for pancake-like grids, 
a 2 = 1, (D) and pencil-like grids a 2 — a\, ( 0 )- The solid line represents the 
theoretically determined values, according to Eq. (3). Error bars are ±cr, where a 
is the standard deviation about the time average. 





264 


A. Scotti tt al 



FIGURE 2. Eddy-viscosity models on pancake-like grids (256 x 16 x 16). (a) kinetic 

energy as a function of time for dynamic model ( ), modified Smagorinsky 

( ) and Smagorinsky- Deardoff ( ). (b) skewness in the worst resolved 

direction, same symbols as in (a). 

perform LES on anisotropic grids characterized by aspect ratios a\ and <22. The 
results are cast in terms of /(ai,a2), by writing 


fdyn( a l 1 ^ 2 ) — 


/ < LijMij > 

0.023 -1 / 2 

y 2 < MijMij > 

A € q 


Figure la shows the time evolution of /dyn(ai,ci2) for three cases: an isotropic grid 
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on 32 3 modes, a pancake-like grid using a 256 x 16 x 16 grid, and a pencil-like grid 
using 128 x 128 x 16 modes. 

In the same way we have computed the time averages of /d yn for aspect ratios 
varying from 1/2 to 1/16. They are plotted in Figure lb together with the value 
obtained from Eq. (3). We see that the dynamic model reproduces the correct trend 
for pancake-like grids, but fails with pencil-like grids. To examine the simulations 
more closely, we now focus on two extreme cases: a 256 x 16 x 16 grid (pancake) 
and a 128 x 128 x 16 grid (pencil). For each case, we compare the dynamic model 
with predictions of the non-dynamic Smagorinsky model and with the non-dynamic 
model but including the correction of Eq. (3). 

4.1 Pancake-like 

Figure 2 shows the total kinetic energy versus time for the three models consid- 
ered. We see that the three models agree quite well. Also, the skewness in the least 
resolved direction does not show marked differences. We conclude that at the large- 
scale level, there is no impact on the model variations even at this high level of grid 
anisotropy. Next, we consider the behavior near the grid scale. The premultiplied 
1-D spectrum is shown in Fig. 3. The traditional Smagorinsky- Deardoff case shows 
a strong peak at wavenumber k i ~ 10. The modified Smagorinsky case remains 
constant at small wavenumbers and dies out at high wavenumbers without showing 
any pile-up. The dynamic model falls somewhere in between, but the value is higher 
than the expected value of C/ v . All models show a rapid decay at wavenumbers 
above 10. 

The fact that all three models decay for k\ > 10 means that those modes that 
cannot have access to all the local triadic interactions experience a high drain of en- 
ergy so that they do not display a Kolmogorov scaling. It appears unlikely that any 
modification of a scalar eddy- viscosity model could compensate for this behavior. 

The analysis of the derivative skewness in the well-resolved direction shows no 
real difference. 

4.2 Pencil-like 

As already mentioned, the dynamic model gives a value for /a yn which is smaller 
than one, in contrast with the theoretical expression, which implies that / must be 
bigger than one. If we look at the large-scale parameters of the flow, energy and 
skewness in the least resolved direction (Fig. 4) we see that the three models again 
give similar answers; note the small value of the skewness in the worst resolved 
direction. But if we consider parameters that are more sensitive to the small scale 
behavior, we notice marked differences. For the dynamic model the Kolmogorov 
constant is too large, about twice as much as expected (Fig. 5). Therefore, the 
“underestimation” of / brings consequences that cannot be ignored at the scales 
near the least resolved direction. Again, scales between the least and best resolved 
directions are much less energetic than the Kolmogorov spectrum, as is clear from 
the rapid drop of the premultiplied spectrum above k\ = 16. On the other hand, the 
modified Smagorinsky model gives too small a value, probably due to overdamped 
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FIGURE 3. Eddy-viscosity models on a pancake-like grid, (a) premultiplied 1-D 

spectrum: dynamic model ( ), modified Smagorinsky ( ) and Smagorinsky- 

Deardoff ( ). (b) derivative skewness in the best resolved direction, same 

symbols as in (a). 

modes near k ~ Finally, the skewness in the best resolved direction is consistent 
with these differences: the smaller the skewness is in magnitude, the more the energy 
piles up. 

4.3 Discussion 

The strongest discrepancy between the theoretically and dynamically determined 
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Figure 4. Eddy- viscosity models on pencil-like grid (128 x 128 x 16). (a) energy 

as a function of time for dynamic model ( ), modified Smagorinsky ( ) 

and Smagorinsky-Deardoff ( ). (b) derivative skewness in the worst resolved 

direction, same symbols as above. 


f(a i,a 2 ) was observed for the case of highly pencil-like grids. For this case, the 
premultiplied spectrum of the dynamic model case showed considerable pile-up, as 
evidenced by much higher values of C(k \ ). In order to understand the causes of this 
behavior, we recall that the dynamic model computes L by sampling the turbulence 
between grid and test filter. It could be argued that for pencil-like grids these modes 
behave essentially as 2D turbulence, with the vorticity aligned in the x 3 direction 
and a concomitant change in the dynamics. To focus on the relevant scales, we 
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t 

FIGURE 5. Pencil-like grid, (a) compensated 1-D spectrum: dynamic model 

( ), modified Smagorinsky ( ) and Smagorinsky-Deardoff ( ). (b) 

derivative skewness in the best resolved direction, same symbols as above. 


have analyzed the vorticity band-pass filtered between test and grid filter (i.e. the 
statistics of = u? — tJ). We find that the variances are not isotropic, and that 
~ u> 2 2 /^ 3 2 ~ 0.75, i.e. the flow is not quite 3-D but not 2-D either. More 
directly related to the small value of L or /d yn obtained from the dynamic model, in 
Fig. 6 we show the PDF of (solid line). The curve is almost symmetrically 

distributed around the origin, and the average value, while positive, is very small 
(< LijMjj >= 4.80). LijMij can be regarded as a measure of energy transfer from 
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large to small scales, with negative values meaning energy backscatter. If we now 
compute the same PDF but using an isotropic test filter at a scale 2A 3 in all three 
directions, we see that the shape of the PDF changes, being now skewed to the 
right (symbols in Fig. 6). The mean value is now < LjjMij >= 31.66. Therefore, 
by sampling larger scales that are more isotropic, the dynamics of the energy transfer 
changes noticeably. 

This observation suggests that in order to improve the performance of the dy- 
namic model in such extreme cases of grid anisotropy, it may be advisable to use a 
test filter which is isotropic, with a length scale twice as large as the worst resolved 
scale. In this case, the grid and test anisotropies differ, and this must be taken 
into account explicitly in the dynamic model formulation. We now implement the 
dynamic model with Eq. (4b) for M tJ , using the expression given in Eq. (3) for 
f(a \ , a 2 ) and /(dj , d 2 )• Using this formulation on a 128 x 128 x 16 simulation yields 
the result shown in Fig. 6. The time trace of / (Fig. 6) shows that it oscillates 
around an average value of 1.44 ± .067, much closer to the expected value of 1.34 
than the value of 0.8 obtained with pencil-like test filtering. At large scales the 
difference between this run and the previous one is small. On the other end, at 
small scales the situation changes as now the premultiplied spectrum (Fig. 7) lies 
flat at 1.4 for k t < 10, very close to the expected value for C K . The skewness in 
the best resolved direction agrees well with the one calculated from the modified 
Smagorinsky model. 

5. Conclusions 

We have run several LES of forced isotropic turbulence on anisotropic grids, us- 
ing three different Smagorinsky models. All three models are able to satisfactorily 
reproduce the very large scales of the flow. This result confirms the general ro- 
bustness of the dynamic model even for the extreme cases considered in this work 
(see Jimenez (1995) for further observations on the dynamic model’s robustness). 
However, none of the models considered is able to give a correct representation 
of the scales smaller than the worst resolved direction, where spectra are strongly 
damped below Kolmogorov values. This is probably due to the fact that the trans- 
fer of energy at very small scales is affected by the lack of similar modes in one 
or more directions. For a related study on the effect of grid anisotropy on velocity 
components and stress anisotropy, see Kaltenbach (1996). 

For the model performance at scales near the cut-off in the worst resolved direc- 
tion, we need to distinguish between pancake grids and pencil grids. For pancake- 
like grids, the non-dynamical Smagorinsky model modified after Scotti et al. (1993) 
and the dynamic model give reasonably good results, while the conventional Smagorin- 
sky model using the Deardoff prescription for A eq shows excessive pile up of energy 
at scales close to the largest mesh size. The anisotropy factor computed from the 
dynamic model shows an increasing trend with anisotropy in accord with the theo- 
retical prediction, although the numerical value is somewhat smaller. For pencil-like 
grids, the Smagorinsky- Deardoff model as well as the modified version give good 
results, with the modified version yielding slightly better results. On the other 


270 


A. Scotti et al 




FIGURE 6. (a) PDF of L tJ Mij computed with same grid but different test filters. 

Both statistics were performed on the same fields simulated on a 128 x 128 x 16 
grid and with test filter cutting off at k{ = l/2k t . The solid line refers to L tj Mij 
computed as in the simulation, while the symbols refer to computed with a 

test filter cutting off at k t = l/2Ar 3 . (b) anisotropy factor / dyn computed with an 

anisotropic test filter ( ) arid with an isotropic (larger scale) test filter ( ). 

The predicted value is 1.34. 


hand, the dynamic model exhibits insufficient dissipation of energy, as shown by 
the fact that the anisotropy factor / dyn becomes smaller than one, and reflected in 
that small scales have excessive energy as compared to the Kolmogorov value. 
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FIGURE 7. (a) compensated 1-D spectrum of dynamic model on pencil-like grid 

with isotropic test filter, (b) derivative skewness in the best resolved direction for 
dynamic model with isotropic test filter ( ) and modified Smagorinsky-model 

( )• 

It would appear that in this particular case the strength of the dynamic model 
becomes its weak point. The dynamic model computes the unknown factor from 
information derived from the smallest resolved scales. But in the case of highly 
anisotropic grids, these scales experience a dynamic which is different from the 
usual one due to the missing modes at large wavenumbers. This in turn affects the 
resolved non-linear interactions embodied in the term , which is what the 

dynamic model samples. Specifically, the number of events during which energy 
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is transferred forward is decreased, which could actually be explained by a partial 
2-dimensionalization of the flow at these scales. 

A proposed improvement is to move the test filter towards larger scales, where 
the combination of more energetic modes and more realistic triadic coupling allows 
a more faithful representation of how energy is exchanged. Indeed, simulations done 
with an isotropic test filter at twice the worst resolved scale show improved results. 
Perhaps not surprisingly, this conclusion is similar to one reached by others in the 
context of dynamic LES using non-spectral numerical methods, such as low-order 
finite differences. There, it has been found advisable to “prefilter” the results and 
shift the test filter to larger scales (Ferziger 1996, Lund 1996) so that the dynamic 
model is not strongly affected by numerical errors occurring near the grid scale. 
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Appendix A 

We assume that the computational domain is covered by A r i points and 

i,j and 1 are unit vectors in the ;r ? y, and z directions. It is well known (see Canuto 
et al (1987)) that the pseudospectral treatment of a 3-D convolution product 
V , , a(mWn) introduces an error. If we denote with w k the true convolution 

Z— /m + n=k \ / v / . . , . 

product and with Wu the calculated one, the following relation holds: 

7 

= w k + 

j= 1 


where the seven extra terms have the form 

Wj = ^2 a(m)b(n) 

m+n = k +ej 


and 

ej = ±Ni i, e 2 = ±iV 2 j, e 3 = ±N 3 1, 
e 4 = ±N\i± N 2 h e 5 = ±Nii± N 3 1, e 6 = ±N 3 1 ± N 2 j, 
e7 = ±N 3 i ± A^ 2 j ± N 3 1. 

The last four terms, (double and triple aliased) can be set to zero if we adopt an 
elliptical truncation, i.e. , if we set to zero all the inodes such that 



The proof is by inspection. 

To remove the single aliased terms we can resort to phase shift. If we premultiply 
all the modes by a factor e lk e , 9 e [0, 2n] x [0, 27 t] x [0, 27r], compute the convolution 
sum and multiply the result by e _ ' k e , the aliased terms now are c ±ie > N >W } , j = 
1,2,3, i.e. we have shifted their phase by an amount ±9jNj. If we do the same 
thing one more time, but this time 9 — > 9+{n /N\ , ir/./V 2 , n/N 3 ) and take the average 
of the results, the aliased terms, being out of phase, will cancel exactly. However, 
this requires doubling the number of FFT’s required for each term to be dealiased. 
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Rogallo (1977) showed that for a multistep scheme such as even-order Runge-Kutta, 
it is possible to control the growth of aliasing essentially at no extra cost. Indeed, 
let us consider the typical step of a 2nd order Runge-Kutta: 


« w+1 = «" + y(*i + F 2 ) 


with F t s being the non-linear terms evaluated recursively. It is important to notice 
that to Oth order in A t they are identical. Therefore, if F\ is evaluated with a shift 9 
and F i with shift 9+(it fN\, Tr/iV^, n/N$), their sum to Oth order is dealiased, leaving 
possibly a contribution to first order. Therefore, the global effect of aliasing is 
pushed to second order. Choosing 9 randomly at each time step further ensures that 
the error does not accumulate over time. Nevertheless the RK-2 method requires 
doubling the FFT’s for each time step. 

In our computation we have used an AB2 scheme, which schematically can be 
written as 


n + l 


u n + ^(3 F n -F n ~ l ) 


with obvious meaning of the symbols. Although to Oth order the alias terms are 
identical in F n and F”" 1 , it is clear that there is no way in which a combination 
of phase shifts can cancel them exactly, since the equation 


Ze ioN - t' 0N = 0 


does not have solutions for a,/? £ [ 0,24 

However, by successive phase-shifts it is still possible to ensure that the error 
does not accumulate. If n is even, the shift is chosen randomly; if n is odd, the shift 
is chosen to be the shift of the previous time step plus (iv/Nx , */N 2l tt/N 3 ). After 
m time steps, the solution can be written as 

u n+m = u n + ^[3(F n + F n+1 + F n+2 -| b F" +ra ) 

- (F"~ 1 -l- F" + F n+1 -| -(- F n+m-1 )]. 

In the two bracketed sums, to the lowest order, all but a few aliased terms (typically 
the first and/or the last) cancel out. This proves that the error does not accumulate, 
and that after m steps the aliasing is still O(At), no matter how big m is. Again, 
the randomness prevents accumulation at higher orders. We have compared results 
obtained with this dealiasing technique with results obtained by zero padding (2-rule 
in the worst resolved direction) without finding any noticeable difference. 



